//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Moment_to_Discrete.cc
 * \author Jeremy Roberts
 * \date   Jul 1, 2011
 * \brief  Moment_to_Discrete member definitions and explicit instantiations.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev::                                               $:Rev of last commit
// $Author::                                            $:Author of last commit
// $Date::                                              $:Date of last commit
//---------------------------------------------------------------------------//

#include <cmath>

#include "angle/Moment_to_Discrete.hh"
#include "angle/Spherical_Harmonics.hh"

#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"

namespace slabtran
{

// Constructor.
template <class D>
Moment_to_Discrete<D>::Moment_to_Discrete(const size_type legendre_order) :
  d_legendre_order(legendre_order),
  d_number_moments(Moments<D>::number_moments(legendre_order))
{
  Require(d_legendre_order >= 0);
  Require(d_number_moments > 0);
}

// Build the operator.
template <class D>
void Moment_to_Discrete<D>::build(SP_angularmesh angularmesh)
{
  Require(angularmesh);

  // Store the angular mesh.
  d_angularmesh = angularmesh;

  // Clear the existing M matrix.
  d_M.clear();
  Check (d_M.empty());

  // Resize M.
  d_number_angles = angularmesh->number_angles();
  d_M.resize(d_number_angles, M_Row(d_number_moments));
  Check (d_M.size() == angularmesh->number_angles());

  // Build the moment-to-discrete operator by looping through each octant
  // and then the angles in each octant.  The ordering is determined by the
  // order of octants in AngularMesh.
  int index = 0;
  // iterators
  typename AngularMesh<D>::const_iterator octant_itr = angularmesh->begin();
  for (; octant_itr != angularmesh->end(); octant_itr++)
  {
      // Get SP to the octant.
      typename AngularMesh<D>::SP_octant octant = *octant_itr;
      Check (octant);
      // Loop over angles in the octant.
      typename Octant<D>::angle_iterator angles = octant->angle_begin();
      for (; angles != octant->angle_end(); angles++)
      {
          // Build M a row at a time, one row for each discrete angle.
          calc_row(index,*angles);
          // Advance the index.
          index++;
      }
  }
  Ensure (angularmesh->number_angles() == index);
}

// 3-d
template <class D>
void Moment_to_Discrete<D>::calc_row(const int angle, const Omega &omega)
{
  Require (angle < d_number_angles);

  // reference to current row
  M_Row &row = d_M[angle];

  // calculate the moments and add them to the row
  for (int l = 0; l <= d_legendre_order; l++)
  {
      double norm = (2.0 * l + 1.0) * Norm<D>::norm;
      for (int m = -l; m <= l; m++)
      {
          // get cardinal moment index
          int i = Moments<D>::index(l,m);
          Check (i < d_number_moments);
          // add harmonic
          row[i] = norm * Spherical_Harmonics::
                   Y_lm(l, m, omega[0], omega[1], omega[2]);
      }
  }
}

// 2-d
template <>
void Moment_to_Discrete<_2D>::calc_row(const int angle, const Omega &omega)
{
  Require (angle < d_number_angles);

  // reference to current row
  M_Row &row = d_M[angle];

  // compute the direction cosine w/r to polar axis
  double xi = sqrt(1.0-omega[0]*omega[0]-omega[1]*omega[1]);

  // loop through l>0 moments and add
  for (int l = 0; l <= d_legendre_order; l++)
  {
    double norm = (2.0 * l + 1.0) * Norm<_2D>::norm;
    for (int m = -l; m <= l; m += 2)
    {
      int i = Moments<_2D>::index(l, m);
      row[i] = norm * Spherical_Harmonics::Y_lm(l, m, omega[0], omega[1], xi);
    }
  }
}

// 1-d
template <>
void Moment_to_Discrete<_1D>::calc_row(const int angle, const Omega &omega)
{
  Require (angle < d_number_angles);

  // reference to current row
  M_Row &row = d_M[angle];

  // calculate the moments and add them to the row
  for (int l = 0; l <= d_legendre_order; l++)
  {
    double norm = (2.0 * l + 1.0) * Norm<_1D>::norm;
    int i = Moments<_1D>::index(l, 0);
    row[i] = norm * Spherical_Harmonics::Y_lm(l, omega);
  }
}

// explicit instantiations
template class Moment_to_Discrete<_1D>;
template class Moment_to_Discrete<_2D>;
template class Moment_to_Discrete<_3D>;


} // end namespace slabtran

//---------------------------------------------------------------------------//
//              end of Moment_to_Discrete.cc
//---------------------------------------------------------------------------//
